## ----------------------------------------------------------------------------
## Title: 04_quantitative.R
## Author: Elsa Voytas
## Last updated: April-13-2025
## ----------------------------------------------------------------------------

library(ggplot2)
library(dplyr)
library(readr)
library(did)
library(gghighlight)
library(did)
library(texreg)
library(tidyverse)
library(lubridate)
library(forcats)
library(stringr)
library(stringi)
library(purrr)
library(tidyr)
library(tibble)
library(patchwork)
library(HonestDiD)
library(staggered)
library(PanelMatch)
library(zoo)
library(bacondecomp)
library(scales)
library(survival)
library(survminer)

source("Functions/reps_functions.R")
panel <- read_csv("Input/panel.df.csv") # This file contains private, personal data and is not included in replication materials
set.seed(5432)
RNGkind("Mersenne-Twister", "Inversion")

#-------------------------------------------------------------------------------
# MAIN PAPER RESULTS
# Producing main quantitative analyses
# Figures 4-8
# ------------------------------------------------------------------------------

## overall ATT
att <- att_gt(yname = "vote",
                   gname = "firstrep",
                   idname = "ID",
                   tname = "time.id",
                   xformla = ~age + female,
                   control_group = "notyettreated",
                   data = panel,
                   est_method = "dr",
                   clustervars = c("ID"))

es <- aggte(att, type = "dynamic", na.rm = T, min_e = -12, max_e = 12)

# Figure 4 ("Estimates of Reparations on Voter Registration for Surviving Victims Approved for Reparations. Coefficient Plot.")
figure4 <- es_plot(es)

ggsave("Output/figure4.pdf", width = 10, height = 5, units="in")
print("Figure 4 complete")

## ATT by age groups
old <- panel %>%
  filter(age.bucket==1)
mid <- panel %>%
  filter(age.bucket==2)
young <- panel %>%
  filter(age.bucket==3)

old.att <- att_gt(yname = "vote",
              tname = "time.id",
              idname = "ID",
              gname = "firstrep",
              xformla = ~age+ female,
              data = old,
              control_group = "notyettreated")

old.es <- aggte(old.att, type = "dynamic", na.rm=T,min_e = -12,max_e = 12)

mid.att <- att_gt(yname = "vote",
              tname = "time.id",
              idname = "ID",
              gname = "firstrep",
              xformla = ~age+ female,
              data = mid,
              control_group = "notyettreated")

mid.es <- aggte(mid.att, type = "dynamic", na.rm=T,min_e = -12,max_e = 12)

young.att <- att_gt(yname = "vote",
                tname = "time.id",
                idname = "ID",
                gname = "firstrep",
                xformla = ~age+ female,
                data = young,
                control_group = "notyettreated")

young.es <- aggte(young.att, type = "dynamic", na.rm=T,min_e = -12,max_e = 12)


old.coefs <- make_coefs_dfs(old.es, "Older")
mid.coefs <- make_coefs_dfs(mid.es, "Middle")
young.coefs <- make_coefs_dfs(young.es, "Younger")

# Figure 5 ("Estimates of Reparations on Voter Registration by Age Group. Coefficient Plot.")
figure5 <- make_coef_plot(old.coefs, mid.coefs, young.coefs)

ggsave("Output/figure5.pdf", width = 12, height = 6, units="in")
print("Figure 5 complete")

## ATT by geography
urban <- panel %>%
  filter(urbandummy==1)
rural <- panel %>%
  filter(urbandummy==0)

urban.att <- att_gt(yname = "vote",
                tname = "time.id",
                idname = "ID",
                gname = "firstrep",
                xformla = ~age+ female,
                data = urban,
                control_group = "notyettreated")

urban.es <- aggte(urban.att, type = "dynamic", na.rm=T,min_e = -12,max_e = 12)

rural.att <- att_gt(yname = "vote",
                tname = "time.id",
                idname = "ID",
                gname = "firstrep",
                xformla = ~age+ female,
                data = rural,
                control_group = "notyettreated")

rural.es <- aggte(rural.att, type = "dynamic", na.rm=T,min_e = -12,max_e = 12)

urban.coefs <- make_coefs_dfs(urban.es, "Urban")
rural.coefs <- make_coefs_dfs(rural.es, "Rural")

# Figure 6 ("Estimates of Reparations on Voter Registration by Geography. Coefficient Plot.")
figure6 <- make_coef_plot(urban.coefs, rural.coefs)

ggsave("Output/figure6.pdf", width = 12, height = 6, units="in")
print("Figure 6 complete")

## ATT by gender
women <- panel %>%
  filter(female==1)
men <- panel %>%
  filter(female==0)

women.att <- att_gt(yname = "vote",
                     tname = "time.id",
                     idname = "ID",
                     gname = "firstrep",
                     xformla = ~age,
                     data = women,
                     control_group =c("notyettreated"))

women.es <- aggte(women.att, type = "dynamic", na.rm=T,min_e = -12,max_e = 12)

men.att <- att_gt(yname = "vote",
                tname = "time.id",
                idname = "ID",
                gname = "firstrep",
                xformla = ~age,
                data = men,
                control_group = "notyettreated")

men.es<- aggte(men.att, type = "dynamic", na.rm=T,min_e = -12,max_e = 12)

men.coefs <- make_coefs_dfs(men.es, "Men")
women.coefs <- make_coefs_dfs(women.es, "Women")

# Figure 7  ("Estimates of Reparations on Voter Registration by Gender. Coefficient Plot.")
figure7 <- make_coef_plot(men.coefs, women.coefs)

ggsave("Output/figure7.pdf", width = 12, height = 6, units="in")
print("Figure 7 complete")

## ATT by income bucket
lowinc <- panel %>%
  filter(bucket==1)

lowerincome.att <- att_gt(yname = "vote",
                      tname = "time.id",
                      idname = "ID",
                      gname = "firstrep",
                      xformla = ~age+ female,
                      data = lowinc,
                      control_group = "notyettreated")

es.lower <- aggte(lowerincome.att, type = "dynamic", na.rm = T, min_e = -12, max_e = 12)

midinc <- panel %>%
  filter(bucket==2)

midincome.att <- att_gt(yname = "vote",
                    tname = "time.id",
                    idname = "ID",
                    gname = "firstrep",
                    xformla = ~age+ female,
                    data = midinc,
                    control_group = "notyettreated")

es.middle <- aggte(midincome.att, type = "dynamic", na.rm = T, min_e = -12, max_e = 12)

highinc <- panel %>%
  filter(bucket==3)

highincome.att <- att_gt(yname = "vote",
                     tname = "time.id",
                     idname = "ID",
                     gname = "firstrep",
                     xformla = ~age + female,
                     data = highinc,
                     control_group = "notyettreated")

es.high <- aggte(highincome.att, type = "dynamic", na.rm=T, min_e=-12, max_e= 12)

low.coefs <- make_coefs_dfs(es.lower, "Low income")
mid.coefs <- make_coefs_dfs(es.middle, "Middle income")
hi.coefs <- make_coefs_dfs(es.high, "High income")

# Figure 8 ("Estimates of Relative Reparations Amount on Voter Registration. Coefficient Plot.")
figure8 <- make_coef_plot(low.coefs, mid.coefs, hi.coefs) + 
  scale_y_continuous(limits = c(-.1,.22), breaks = seq(-.1, .22, by = 0.02))

ggsave("Output/figure8.pdf", width = 12, height = 6, units="in")
print("Figure 8 complete")

#-------------------------------------------------------------------------------
# APPENDIX RESULTS
# Producing Tables A1-A5, Figures A2-A3
# ------------------------------------------------------------------------------

# Reading in data
full.df <- read_csv("Input/full.df.csv")
voting <- read_csv("Input/voting.csv")
valech <- read_csv("Input/valech.csv")
first.df <- read_csv("Input/sample.df.csv")
panel <- read_csv("Input/panel.df.csv")

# Table A1: Monthly reparations amounts
# Reparations amounts from Ley 19.992. 
# Minimum wage from Ley 20.039. 
# Conversions using the Minneapolis Fed inflation calculator and exchangerates.org.

# Figure A1 ("Trends in reparations approvals")
figurea1 <- full.df %>%
  filter(date_determined< "2010-09-01") %>%
  ggplot(aes(x=date_determined)) +
  geom_histogram(bins=150) +
  theme_minimal() +
  geom_vline(xintercept = as.numeric(as.Date("2005-12-11")), linetype="dashed", 
             color = "red") + #senate and president
  geom_vline(xintercept = as.numeric(as.Date("2006-01-15")), linetype="dashed", 
             color = "blue") + #president run off
  geom_vline(xintercept = as.numeric(as.Date("2008-10-26")), linetype="dashed", 
             color = "red") + #local election
  geom_vline(xintercept = as.numeric(as.Date("2009-12-13")), linetype="dashed", 
             color = "red") + # presidential election
  scale_x_date(date_breaks = "6 months", labels = date_format("%m-%Y")) +
  labs(x="Date of issue", y="Count")

ggsave("Output/figurea1.pdf", width = 10, height = 6, units="in")
print("Figure A1 complete")

# Table A2 ("Sample population")

# First row
top.left <- full.df %>%
  filter(finscrip < min(date_determined)) %>%
  filter(date_determined < "2010-09-01")
nrow(top.left)

top.right <- full.df %>%
  filter(finscrip > min(date_determined) | is.na(finscrip)) %>%
  filter(date_determined < "2010-09-01")
nrow(top.right)

# Second row
mid.left <- full.df %>%
  filter(finscrip < min(date_determined)) %>%
  filter(date_determined > "2010-09-01")
nrow(mid.left)

mid.right <- full.df %>%
  filter(finscrip > min(date_determined) | is.na(finscrip)) %>%
  filter(date_determined > "2010-09-01")
nrow(mid.right)

# Last row
# Reading in and cleaning data on surviving victims who did not claim reparations
no <- read_csv("Input/norepsIPS.csv")
no <- no %>%
  na.omit() # removing empty rows

# adjusting name format for merge with rest of data
not <- no %>% 
  tidyr::separate(Nombre_Nomina_Completo,into=c("lastname1", "lastname2",
                                                "firstname1","firstname2","firstname3",
                                                "firstname4","firstname5","firstname6",
                                                "firstname7"), 
                  sep='\\s')
not$firstname1 <-toupper(not$firstname1)
not$firstname1<-stri_trans_general(not$firstname1,"Latin-ASCII")
not$firstname2 <-toupper(not$firstname2)
not$firstname2<-stri_trans_general(not$firstname2,"Latin-ASCII")
not$firstname3 <-toupper(not$firstname3)
not$firstname3<-stri_trans_general(not$firstname3,"Latin-ASCII")
not$firstname4 <-toupper(not$firstname4)
not$firstname4<-stri_trans_general(not$firstname4,"Latin-ASCII")
not$firstname5 <-toupper(not$firstname5)
not$firstname5<-stri_trans_general(not$firstname5,"Latin-ASCII")
not$firstname6 <-toupper(not$firstname6)
not$firstname6<-stri_trans_general(not$firstname6,"Latin-ASCII")
not$firstname7 <-toupper(not$firstname7)
not$firstname7 <-stri_trans_general(not$firstname7,"Latin-ASCII")
not$lastname1 <-toupper(not$lastname1)
not$lastname1<-stri_trans_general(not$lastname1,"Latin-ASCII")
not$lastname2 <-toupper(not$lastname2)
not$lastname2<-stri_trans_general(not$lastname2,"Latin-ASCII")
not <- not %>% 
  unite(., col = "name",c("lastname1","lastname2", "firstname1", "firstname2",
                          "firstname3","firstname4","firstname5", "firstname6", "firstname7"),
        na.rm=TRUE, sep = " ") %>%
  mutate(name = gsub("\\,","",name))

not <- not %>%
  distinct()

# checking to see if in IPS
not <- subset(not, !(name %in% full.df$name))

# merging with voting data
not.voting <- left_join(not, voting, by=c("name"))
not.voting$vote <- ifelse(!is.na(not.voting$finscrip),1,0)
not.voting$finscrip <- as.character(not.voting$finscrip)
not.voting$finscrip <- as.Date(not.voting$finscrip, format="%Y%m%d")

not.voting$prereps <- ifelse(not.voting$finscrip>min(panel$date_determined)|
                               not.voting$vote==0,0,
                             1)
table(not.voting$prereps)
#  0   1 
# 698 524 

print("Table A2 complete")

## Survival analysis
surv.reps <- first.df %>%
  filter(date_determined<"2010-09-01"&finscrip>"2005-01-04"|
           date_determined<"2010-09-01"&is.na(finscrip)) %>%
  dplyr::select(finscrip) %>%
  mutate(reps=1)

surv.noreps <- not.voting %>%
  filter(finscrip >"2005-01-04"|is.na(finscrip)) %>%
  dplyr::select(finscrip) %>%
  mutate(reps=0)

surv <- rbind(surv.reps, surv.noreps)

surv$time <- surv$finscrip - min(panel$date_determined)
surv$status <- ifelse(is.na(surv$finscrip),0,1)
surv$time <- ifelse(is.na(surv$time), 
                    max(first.df$finscrip,na.rm=T)- min(panel$date_determined),
                    surv$time)

s <- Surv(surv$time, surv$status)
sfit <- survfit(Surv(surv$time,surv$status)~reps, data=surv)

# Figure A2 ("Survival plot of victims receiving versus not receiving reparations payments")
ggsurv <- ggsurvplot(sfit,fun="event",conf.int = T,pval=T,palette="grey",legend="none")+
  labs(x="Time in days",y="Probability of registering to vote")

figurea2 <- ggsurv$plot <- ggsurv$plot +
  ggplot2::annotate(geom="text", x=1350, y=.36, label="Victims with reparations (n=5,432)",
                    color="black")+
  ggplot2::annotate( geom="text", x=1580, y=.18, label="Victims without reparations (n=698)",
                     color="black") +
  theme(axis.line = element_line(colour = "black"),
        panel.grid.major = element_line(colour = "grey90"),
        panel.grid.minor = element_line(colour = "grey90"),
        panel.border = element_blank(),
        panel.background = element_blank())

ggsave("Output/figurea2.pdf", width = 12, height = 6, units="in")
print("Figure A2 complete")

# Table A3: Population and sample characteristics
# Population column data from 2002 Census (INE)
# Age and gender breakdown from: https://redatam-ine.ine.cl/redbin/RpWebEngine.exe/Portal?BASE=CENSO_2002&lang=esp
# Urban/rural breakdown from: https://www.ine.gob.cl/docs/default-source/censo-de-poblacion-y-vivienda/publicaciones-y-anuarios/2002/sintesiscensal-2002.pdf
# Both accessed May 13, 2024 and saved in replication package

# Column 3 -- Voting population
voting$age = 2010 - as.numeric(voting$birthyear)

voting <- voting %>%
  mutate(age_group = case_when(age >=20 & age <= 24~"group2024",
                               age >=25 & age <= 29~ "group2529",
                               age >=30 & age <= 34~ "group3034",
                               age >=35 & age <= 39~ "group3539",
                               age >=40 & age <= 44~ "group4044",
                               age >=45 & age <= 49~ "group4549",
                               age >=50 & age <= 54~ "group5054",
                               age >=55 & age <= 59~ "group5559",
                               age >=60 & age <= 64~ "group6064",
                               age >=65 & age <= 69~ "group6569",
                               age >= 70 & age <= 74 ~ "group7074",
                               age >= 75 & age <= 79 ~ "group7579",
                               age >= 80 ~ "group80+",
                               TRUE~NA_character_))
age_sum <- voting %>%
  group_by(age_group) %>%
  dplyr::summarize(count = n()) %>%
  ungroup() %>%
  mutate(prop = count/sum(count)*100)

age_sum <- age_sum %>%
  mutate(across(where(is.numeric), .fns = function(x) {format(round(x, 2), nsmall = 2)}))

gender_sum <-  voting %>%
  group_by(sexo) %>%
  dplyr::summarize(count = n()) %>%
  ungroup() %>%
  mutate(prop = (count/sum(count)*100))
gender_sum <- gender_sum %>%
  mutate(across(where(is.numeric), .fns = function(x) {format(round(x, 2), nsmall = 2)}))

# Column 4 -- Registered surviving victims
valech.merge <- left_join(valech, voting, by=c("rut","dv"))
valech.merge$finscrip <- as.character(valech.merge$finscrip)
valech.merge$finscrip <- as.Date(valech.merge$finscrip, format="%Y%m%d")
valech.merge$yrinscrip <- format(as.Date(valech.merge$finscrip, format="%Y-%m-%d"),"%Y")

voting_victims <- valech.merge %>%
  filter(!is.na(finscrip)) %>%
  mutate(birthyear = substr(fnacim, 1, 4))
voting_victims$age = 2010 - as.numeric(voting_victims$birthyear)

voting_victims <- voting_victims %>%
  mutate(age_group = case_when(age >=20 & age <= 24~"group2024",
                               age >=25 & age <= 29~ "group2529",
                               age >=30 & age <= 34~ "group3034",
                               age >=35 & age <= 39~ "group3539",
                               age >=40 & age <= 44~ "group4044",
                               age >=45 & age <= 49~ "group4549",
                               age >=50 & age <= 54~ "group5054",
                               age >=55 & age <= 59~ "group5559",
                               age >=60 & age <= 64~ "group6064",
                               age >=65 & age <= 69~ "group6569",
                               age >= 70 & age <= 74 ~ "group7074",
                               age >= 75 & age <= 79 ~ "group7579",
                               age >= 80 ~ "group80+",
                               TRUE~NA_character_))

age_sum <- voting_victims %>%
  group_by(age_group) %>%
  dplyr::summarize(count = n()) %>%
  ungroup() %>%
  mutate(prop = (count/sum(count)*100))
age_sum <- age_sum %>%
  mutate(across(where(is.numeric), .fns = function(x) {format(round(x, 2), nsmall = 2)}))

gender_sum <-  voting_victims %>%
  group_by(sexo) %>%
  dplyr::summarize(count = n()) %>%
  ungroup() %>%
  mutate(prop = (count/sum(count)*100))
gender_sum <- gender_sum %>%
  mutate(across(where(is.numeric), .fns = function(x) {format(round(x, 2), nsmall = 2)}))

# Column 5 -- Reparations recipients
# all reparations recipients
full.df$age.table = 2010 - as.numeric(full.df$birthyr)

full.df <- full.df %>%
  mutate(age_group = case_when(age.table >=20 & age.table <= 24~"group2024",
                               age.table >=25 & age.table <= 29~ "group2529",
                               age.table >=30 & age.table <= 34~ "group3034",
                               age.table >=35 & age.table <= 39~ "group3539",
                               age.table >=40 & age.table <= 44~ "group4044",
                               age.table >=45 & age.table <= 49~ "group4549",
                               age.table >=50 & age.table <= 54~ "group5054",
                               age.table >=55 & age.table <= 59~ "group5559",
                               age.table >=60 & age.table <= 64~ "group6064",
                               age.table >=65 & age.table <= 69~ "group6569",
                               age.table >= 70 & age.table <= 74 ~ "group7074",
                               age.table >= 75 & age.table <= 79 ~ "group7579",
                               age.table >= 80 ~ "group80+",
                               TRUE~NA_character_))

age_sum <- full.df %>%
  group_by(age_group) %>%
  dplyr::summarize(count = n()) %>%
  ungroup() %>%
  mutate(prop = (count/sum(count)*100))
age_sum <- age_sum %>%
  mutate(across(where(is.numeric), .fns = function(x) {format(round(x, 2), nsmall = 2)}))

gender_sum <-  full.df %>%
  group_by(female) %>%
  dplyr::summarize(count = n()) %>%
  ungroup() %>%
  mutate(prop = (count/sum(count)*100))
gender_sum <- gender_sum %>%
  mutate(across(where(is.numeric), .fns = function(x) {format(round(x, 2), nsmall = 2)}))

urban_sum <-  full.df %>%
  group_by(urbandummy) %>%
  dplyr::summarize(count = n()) %>%
  ungroup() %>%
  mutate(prop = (count/sum(count)*100))
urban_sum <- urban_sum %>%
  mutate(across(where(is.numeric), .fns = function(x) {format(round(x, 2), nsmall = 2)}))

# sample
first.df$age.table = 2010 - as.numeric(first.df$birthyr)

first.df <- first.df %>%
  mutate(age_group = case_when(age.table >=20 & age.table <= 24~"group2024",
                               age.table >=25 & age.table <= 29~ "group2529",
                               age.table >=30 & age.table <= 34~ "group3034",
                               age.table >=35 & age.table <= 39~ "group3539",
                               age.table >=40 & age.table <= 44~ "group4044",
                               age.table >=45 & age.table <= 49~ "group4549",
                               age.table >=50 & age.table <= 54~ "group5054",
                               age.table >=55 & age.table <= 59~ "group5559",
                               age.table >=60 & age.table <= 64~ "group6064",
                               age.table >=65 & age.table <= 69~ "group6569",
                               age.table >= 70 & age.table <= 74 ~ "group7074",
                               age.table >= 75 & age.table <= 79 ~ "group7579",
                               age.table >= 80 ~ "group80+",
                               TRUE~NA_character_))

age_sum <- first.df %>%
  group_by(age_group) %>%
  dplyr::summarize(count = n()) %>%
  ungroup() %>%
  mutate(prop = (count/sum(count)*100))
age_sum <- age_sum %>%
  mutate(across(where(is.numeric), .fns = function(x) {format(round(x, 2), nsmall = 2)}))

gender_sum <-  first.df %>%
  group_by(female) %>%
  dplyr::summarize(count = n()) %>%
  ungroup() %>%
  mutate(prop = (count/sum(count)*100))
gender_sum <- gender_sum %>%
  mutate(across(where(is.numeric), .fns = function(x) {format(round(x, 2), nsmall = 2)}))

urban_sum <-  first.df %>%
  group_by(urbandummy) %>%
  dplyr::summarize(count = n()) %>%
  ungroup() %>%
  mutate(prop = (count/sum(count)*100))
urban_sum <- urban_sum %>%
  mutate(across(where(is.numeric), .fns = function(x) {format(round(x, 2), nsmall = 2)}))

print("Table A3 complete")

# Table A4: Profile of surviving victims
# Data from Valech Report, available at: https://www.indh.cl/bb/wp-content/uploads/2017/01/informe.pdf
# Accessed May 13, 2024

## Figure A3: Surviving victim registration trends alongside key developments
# See code in 02_registration.R

## Table A5: Covariate balance
first.df <- first.df %>%
  mutate(before=ifelse(date_determined<median(date_determined),1,0))
options(scipen = 999)
t.test(first.df$female~first.df$before)
t.test(as.integer(first.df$age)~first.df$before)
t.test(as.numeric(first.df$urban)~first.df$before)

print("Table A5 complete")

#-------------------------------------------------------------------------------
# Producing Tables A6-A11 and Figures A4, A17
# ------------------------------------------------------------------------------

# Aggregate Results
simple <- aggte(att, type = "simple", na.rm=T)


## overall ATT -- without covariates 
att2 <- att_gt(yname = "vote",
               gname = "firstrep",
               idname = "ID",
               tname = "time.id",
               #xformla = ~age + female,
               control_group = "notyettreated",
               data = panel,
               est_method = "dr",
               clustervars = c("ID"))

# Producing Table A6: Modeling the Effects of Reparations on Voter Registration (Aggregate)
simple2 <- aggte(att2, type = "simple", na.rm=T)
tr.c <- createTexreg(coef = simple$overall.att, coef.names = "ATT", 
                     se = simple$overall.se)
tr <- createTexreg(coef = simple2$overall.att, coef.names = "ATT", 
                   se = simple2$overall.se)
tablea6 <- texreg(c(tr.c, tr),  digits=3, custom.note= "", 
                  custom.gof.rows = list("95% CIs" = c("0.008,0.129", "0.034,0.171"),
                                         "Covariates" = c("Yes", "No"),
                                         N = c(simple$DIDparams$n, simple2$DIDparams$n)))

tablea6_models <- list(
  model1_with_covariates = list(
    att = simple$overall.att,
    se = simple$overall.se,
    n = simple$DIDparams$n
  ),
  model2_no_covariates = list(
    att = simple2$overall.att,
    se = simple2$overall.se,
    n = simple2$DIDparams$n
  )
)

saveRDS(tablea6_models, "Output/tablea6.rds")

print("Table A6 complete")

# Full Model Results
# Varying baseline
es <- aggte(att, type = "dynamic", na.rm = T, min_e = -12, max_e = 12)

es2 <- aggte(att2, type = "dynamic", na.rm = T, min_e = -12, max_e = 12)

# Universal baseline
att_u <- att_gt(yname = "vote",
                gname = "firstrep",
                idname = "ID",
                tname = "time.id",
                xformla = ~age + female,
                control_group = "notyettreated",
                data = panel,
                base_period = "universal",
                est_method = "dr",
                clustervars = c("ID"))

# without covariates
att_u2 <- att_gt(yname = "vote",
                 gname = "firstrep",
                 idname = "ID",
                 tname = "time.id",
                 #xformla = ~age + female,
                 control_group = "notyettreated",
                 data = panel,
                 base_period = "universal",
                 est_method = "dr",
                 clustervars = c("ID"))

es_u <- aggte(att_u, type = "dynamic", na.rm = T, min_e = -12, max_e = 12)

es_u2 <- aggte(att_u2, type = "dynamic", na.rm = T, min_e = -12, max_e = 12)


# Table A7: Modeling the Effects of Reparations on Voter Registration
col1 <- extractFromDataFrame(tidy(es))
col2 <- extractFromDataFrame(tidy(es2))
col3 <- extractFromDataFrame(tidy(es_u))
col4 <- extractFromDataFrame(tidy(es_u2))
tablea7 <- texreg(c(col1, col2, col3, col4),  digits=3,
                  custom.note= "", caption = "", label = "",
                  custom.gof.rows = list("Covariates" = c("Yes", "No", "Yes", "No"),
                                         Individuals = c(es$DIDparams$n, es2$DIDparams$n,
                                                         es_u$DIDparams$n, es_u2$DIDparams$n),
                                         Baseline = c("Varying", "Varying", "Universal", "Universal")))

tablea7_models <- list(
  model1_varying_covs = list(
    e = es$egt,
    att = es$att.egt,
    se = es$se.egt,
    n = es$DIDparams$n
  ),
  model2_varying_nocovs = list(
    e = es2$egt,
    att = es2$att.egt,
    se = es2$se.egt,
    n = es2$DIDparams$n
  ),
  model3_universal_covs = list(
    e = es_u$egt,
    att = es_u$att.egt,
    se = es_u$se.egt,
    n = es_u$DIDparams$n
  ),
  model4_universal_nocovs = list(
    e = es_u2$egt,
    att = es_u2$att.egt,
    se = es_u2$se.egt,
    n = es_u2$DIDparams$n
  )
)

saveRDS(tablea7_models, "Output/tablea7.rds")

print("Table A7 complete")

# Results by age and geography
# Table A8: Effects of Reparations on Voter Registration by Municipality Income Level
col1 <- extractFromDataFrame(tidy(es.lower))
col2 <- extractFromDataFrame(tidy(es.middle))
col3 <- extractFromDataFrame(tidy(es.high))
tablea8 <- texreg(c(col1, col2, col3),  digits=3,
                  custom.model.names = c("Low", "Middle", "High"),
                  custom.note= "", caption = "", label = "",
                  custom.header = list("Income level" = 1:3),
                  custom.gof.rows = list("Covariates" = c("Yes", "Yes", "Yes"),
                                         Individuals = c(es.lower$DIDparams$n, es.middle$DIDparams$n,
                                                         es.high$DIDparams$n),
                                         Baseline = c("Varying", "Varying", "Varying")))

tablea8_models <- list(
  low_income = list(
    e = es.lower$egt,
    att = es.lower$att.egt,
    se = es.lower$se.egt,
    n = es.lower$DIDparams$n
  ),
  middle_income = list(
    e = es.middle$egt,
    att = es.middle$att.egt,
    se = es.middle$se.egt,
    n = es.middle$DIDparams$n
  ),
  high_income = list(
    e = es.high$egt,
    att = es.high$att.egt,
    se = es.high$se.egt,
    n = es.high$DIDparams$n
  )
)

saveRDS(tablea8_models, "Output/tablea8.rds")

print("Table A8 complete")

# Table A9: Effects of Reparations on Voter Registration by Age
col1 <- extractFromDataFrame(tidy(young.es))
col2 <- extractFromDataFrame(tidy(mid.es))
col3 <- extractFromDataFrame(tidy(old.es))
tablea9 <- texreg(c(col1, col2, col3),  digits=3,
                  custom.model.names = c("Younger", "Middle", "Older"),
                  custom.note= "", caption = "", label = "",
                  custom.header = list("Age bracket" = 1:3),
                  custom.gof.rows = list("Covariates" = c("Yes", "Yes", "Yes"),
                                         Individuals = c(young.es$DIDparams$n, mid.es$DIDparams$n,
                                                         old.es$DIDparams$n),
                                         Baseline = c("Varying", "Varying", "Varying")))

tablea9_models <- list(
  younger = list(
    e = young.es$egt,
    att = young.es$att.egt,
    se = young.es$se.egt,
    n = young.es$DIDparams$n
  ),
  middle = list(
    e = mid.es$egt,
    att = mid.es$att.egt,
    se = mid.es$se.egt,
    n = mid.es$DIDparams$n
  ),
  older = list(
    e = old.es$egt,
    att = old.es$att.egt,
    se = old.es$se.egt,
    n = old.es$DIDparams$n
  )
)

saveRDS(tablea9_models, file = "Output/tablea9.rds")

print("Table A9 complete")

# Table A10: Effects of Reparations on Voter Registration by Geography

col1 <- extractFromDataFrame(tidy(urban.es))
col2 <- extractFromDataFrame(tidy(rural.es))
tablea10 <- texreg(c(col1, col2),  digits=3,
                   custom.model.names = c("Urban", "Rural"),
                   custom.note= "", caption = "", label = "",
                   custom.gof.rows = list("Covariates" = c("Yes", "Yes"),
                                          Individuals = c(urban.es$DIDparams$n, rural.es$DIDparams$n),
                                          Baseline = c("Varying", "Varying")))

tablea10_models <- list(
  urban = list(
    e = urban.es$egt,
    att = urban.es$att.egt,
    se = urban.es$se.egt,
    n = urban.es$DIDparams$n
  ),
  rural = list(
    e = rural.es$egt,
    att = rural.es$att.egt,
    se = rural.es$se.egt,
    n = rural.es$DIDparams$n
  )
)

# Save to RDS
saveRDS(tablea10_models, file = "Output/tablea10.rds")

print("Table A10 complete")

# Table A11: Effects of Reparations on Voter Registration by Gender
col1 <- extractFromDataFrame(tidy(men.es))
col2 <- extractFromDataFrame(tidy(women.es))
tablea11 <- texreg(c(col1, col2),  digits=3,
                   custom.model.names = c("Men", "Women"),
                   custom.note= "", caption = "", label = "",
                   custom.gof.rows = list("Covariates" = c("Yes", "Yes"),
                                          Individuals = c(men.es$DIDparams$n, women.es$DIDparams$n),
                                          Baseline = c("Varying", "Varying")))
tablea11_models <- list(
  men = list(
    e = men.es$e,
    att = men.es$att,
    se = men.es$se,
    n = men.es$DIDparams$n
  ),
  women = list(
    e = women.es$e,
    att = women.es$att,
    se = women.es$se,
    n = women.es$DIDparams$n
  )
)

# Save to RDS file
saveRDS(tablea11_models, file = "Output/tablea11.rds")

print("Table A11 complete")

## Figure A4 ("Estimates of Reparations on Voter Registration for Surviving
# Victims Approved for Reparations. Universal Baseline. Coefficient Plot.")
figurea4 <- es_plot(es_u)

ggsave("Output/figurea4.pdf", width = 10, height = 5, units="in")
print("Figure A4 complete")

## Figure A5 ("Estimates of Reparations on Voter Registration for Surviving
# Victims Approved for Reparations. Balanced Sample. Coefficient Plot.")
es.balanced <- aggte(att, type = "dynamic", na.rm = T, balance_e = 12)
figurea5 <- es_plot(es.balanced)

ggsave("Output/figurea5.pdf", width = 10, height = 5, units="in")
print("Figure A5 complete")

## Figure A6 ("Estimates of Reparations on Voter Registration for Surviving
# Victims Approved for Reparations. Never-treated Control Group. Coefficient Plot.")

# Never-treated group - this includes people who claimed reparations AFTER my voting data ends
# these folks now serve as a "never-treated" counterfactual
never_sample <- full.df[which(is.na(full.df$finscrip)|full.df$finscrip>min(full.df$date_determined)),]
never_sample <- never_sample %>%
  mutate(EndDate = as.Date(max(finscrip,na.rm=T)),
         StartDate = as.yearmon(StartDate, format = "%Y-%m"),
         EndDate = as.yearmon(EndDate, format = "%Y-%m"))

panel_n <- never_sample %>%
  dplyr::group_by(ID, F_Nac) %>%
  mutate(time = list(seq.Date(as.Date(StartDate), as.Date(EndDate), by = "month"))) %>%
  tidyr::unnest()

# creating treatment and outcome variables
panel_n <- panel_n %>%
  mutate(vote = ifelse(finscrip>time,0,1)) # if voter registration date is after row time, code 0, otw 1
panel_n$vote[is.na(panel_n$vote)] <- 0 # if NA, that means registration date is NA (unregistred), so code as 0
panel_n$rep.date <- panel_n$date_determined 
panel_n <- panel_n %>%
  mutate(reparation = ifelse(rep.date > time, 0,1))

# create time ID (numeric variable for use in packages)
times <- data.frame(time = sort(unique(panel_n$time)),
                    time.id = 1:length(unique(panel_n$time)))
panel_n.df <- merge(panel_n, times, by = "time") %>%
  distinct(name, time.id, .keep_all = T) %>%
  group_by(name) %>%
  mutate(firstrep = time.id[which.max(reparation)]) # variable
panel_n.df$firstrep <- ifelse(panel_n.df$date_determined>"2010-09-01",Inf,
                              panel_n.df$firstrep) # convert afters to Inf


att_n_1 <- att_gt(yname = "vote",
                  tname = "time.id",
                  idname = "ID",
                  gname = "firstrep",
                  xformla = ~age + female,
                  data = panel_n.df,
                  control_group = c("nevertreated"))

es_n_1 <- aggte(att_n_1, type = "dynamic", na.rm = T, min_e = -12, max_e = 12)
figurea6 <- es_plot(es_n_1)

ggsave("Output/figurea6.pdf", width = 10, height = 5, units="in")
print("Figure A6 complete")

## Figure A7 ("Estimates of Reparations on Voter Registration for Surviving
# Victims Approved for Reparations. Longer pre-treatment period. Coefficient Plot")
es.prelong <- aggte(att_u, type = "dynamic", na.rm = T, min_e = -48, max_e = 12)
figurea7 <- es_plot_pre(es.prelong) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

ggsave("Output/figurea7.pdf", width = 10, height = 5, units="in")
print("Figure A7 complete")

## Figure A8 ("Estimates of Reparations on Voter Registration for Surviving Victims Approved for
# Reparations. Longer post-treatment period. Coefficient Plot.")
es.postlong <- aggte(att, type = "dynamic", na.rm = T, min_e = -12, max_e = 48)
figurea8 <- es_plot_post(es.postlong) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

ggsave("Output/figurea8.pdf", width = 10, height = 5, units="in")
print("Figure A8 complete")

## Figure A9 ("Estimates of Reparations on Voter Registration for Surviving Victims Aprpoved for
# Reparations. First year only. Coefficient Plot.")
panel.yr1 <- panel[which(panel$firstrep<=12),]

att.yr1 <- att_gt(yname = "vote",
                  gname = "firstrep",
                  idname = "ID",
                  tname = "time.id",
                  xformla = ~age + female,
                  control_group = "notyettreated",
                  data = panel.yr1,
                  est_method = "dr",
                  clustervars = c("ID"))

es.yr1 <- aggte(att, type = "dynamic", na.rm = T, min_e = -12, max_e = 12)
figurea9 <- es_plot(es.yr1)

ggsave("Output/figurea9.pdf", width = 10, height = 5, units="in")
print("Figure A9 complete")

## Figure A10 ("Estimates of Reparations on Voter Registration for Surviving Victims Approved for
# Reparations. Calendar time. Coefficient Plot.")
calendar.es <- aggte(att, na.rm=T,type = "calendar")
figurea10 <- ggdid(calendar.es)

ggsave("Output/figurea10.pdf", width = 16, height = 8, units="in")
print("Figure A10 complete")

# ------------------------------------------------------------------------------
# Sensitivity analysis
# ------------------------------------------------------------------------------

## 12 pre-treatment periods
hd_1 <- honest_did(es = es_u, 
                   e = 1,
                   type="relative_magnitude", grid.ub = .6, grid.lb=-.6)

hd_1$robust_ci <- hd_1$robust_ci[-1,]
hd1plot <- createSensitivityPlot_relativeMagnitudes(hd_1$robust_ci,
                                                    hd_1$orig_ci) +
  ggtitle("e=1") + theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_2 <- honest_did(es = es_u, 
                   e = 2,
                   type="relative_magnitude", grid.ub = .6, grid.lb=-.6)
hd_2$robust_ci <- hd_2$robust_ci[-1,]
hd2plot <- createSensitivityPlot_relativeMagnitudes(hd_2$robust_ci,
                                                    hd_2$orig_ci) +
  ggtitle("e=2") + theme_minimal() + theme(legend.position="bottom") + scale_color_discrete(name="")
hd_3 <- honest_did(es = es_u, 
                   e = 3,
                   type="relative_magnitude",grid.ub = .6, grid.lb=-.6)
hd_3$robust_ci <- hd_3$robust_ci[-1,]
hd3plot <- createSensitivityPlot_relativeMagnitudes(hd_3$robust_ci,
                                                    hd_3$orig_ci) +
  ggtitle("e=3") + theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")
hd_4 <- honest_did(es = es_u, 
                   e = 4,
                   type="relative_magnitude", grid.ub = .6, grid.lb=-.6)
hd_4$robust_ci <- hd_4$robust_ci[-1,]
hd4plot <- createSensitivityPlot_relativeMagnitudes(hd_4$robust_ci,
                                                    hd_4$orig_ci) +
  ggtitle("e=4") + theme_minimal() + theme(legend.position="bottom") + scale_color_discrete(name="")
hd_5 <- honest_did(es = es_u, 
                   e = 5,
                   type="relative_magnitude", grid.ub = .6, grid.lb=-.6)
hd_5$robust_ci <- hd_5$robust_ci[-1,]
hd5plot <- createSensitivityPlot_relativeMagnitudes(hd_5$robust_ci,
                                                    hd_5$orig_ci) +
  ggtitle("e=5") + theme_minimal() + theme(legend.position="bottom") + scale_color_discrete(name="")
hd_6 <- honest_did(es = es_u, 
                   e = 6,
                   type="relative_magnitude", grid.ub = .6, grid.lb=-.6)
hd_6$robust_ci <- hd_6$robust_ci[-1,]
hd6plot <- createSensitivityPlot_relativeMagnitudes(hd_6$robust_ci,
                                                    hd_6$orig_ci)+
  ggtitle("e=6")+  theme_minimal()+theme(legend.position="bottom") + scale_color_discrete(name="")
hd_7 <- honest_did(es = es_u, 
                   e = 7,
                   type="relative_magnitude", grid.ub = 3, grid.lb=-1.8)
hd_7$robust_ci <- hd_7$robust_ci[-1,]
hd7plot <- createSensitivityPlot_relativeMagnitudes(hd_7$robust_ci,
                                                    hd_7$orig_ci) +
  ggtitle("e=7")+ theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="") 

hd_8 <- honest_did(es = es_u, 
                   e = 8,
                   type="relative_magnitude", grid.ub = 3, grid.lb=-1.8)
hd_8$robust_ci <- hd_8$robust_ci[-1,]
hd8plot <- createSensitivityPlot_relativeMagnitudes(hd_8$robust_ci,
                                                    hd_8$orig_ci) +
  ggtitle("e=8")+ theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")
hd_9 <- honest_did(es = es_u, 
                   e = 9,
                   type="relative_magnitude", grid.ub = 3, grid.lb=-1.8)
hd_9$robust_ci <- hd_9$robust_ci[-1,]
hd9plot <- createSensitivityPlot_relativeMagnitudes(hd_9$robust_ci,
                                                    hd_9$orig_ci) +
  ggtitle("e=9")+ theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")
hd_10 <- honest_did(es = es_u, 
                    e = 10,
                    type="relative_magnitude", grid.ub = 3, grid.lb=-1.8)
hd_10$robust_ci <- hd_10$robust_ci[-1,]
hd10plot <- createSensitivityPlot_relativeMagnitudes(hd_10$robust_ci,
                                                     hd_10$orig_ci) +
  ggtitle("e=10")+ theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")
hd_11 <- honest_did(es = es_u, 
                    e = 11,
                    type="relative_magnitude", grid.ub = 3, grid.lb=-1.8)
hd_11$robust_ci <- hd_11$robust_ci[-1,]
hd11plot <- createSensitivityPlot_relativeMagnitudes(hd_11$robust_ci,
                                                     hd_11$orig_ci) +
  ggtitle("e=11")+ theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")
hd_12 <- honest_did(es = es_u, 
                    e = 12,
                    type="relative_magnitude", grid.ub = 3, grid.lb=-1.8)
hd_12$robust_ci <- hd_12$robust_ci[-1,]
hd12plot <- createSensitivityPlot_relativeMagnitudes(hd_12$robust_ci,
                                                     hd_12$orig_ci) + 
  ggtitle("e=12")+ theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

combined <- (hd1plot + hd2plot + hd3plot) /
  (hd4plot + hd5plot + hd6plot)/  (hd7plot +
                                     hd8plot + hd9plot) / (hd10plot+ hd11plot +
                                                             hd12plot)  & theme(legend.position = "bottom")
figurea11 <- combined + plot_layout(guides = "collect")

ggsave("Output/figurea11.pdf", width = 8, height = 8, units="in")
print("Figure A11 complete")

## 4 pre-treatment periods

es_u4 <- aggte(att_u, type = "dynamic", na.rm = T, min_e = -4, max_e = 12)
hd_1 <- honest_did(es = es_u4, 
                   e = 1,
                   type="relative_magnitude")
hd_1$robust_ci <- hd_1$robust_ci[-1,]

hd1plot <- createSensitivityPlot_relativeMagnitudes(hd_1$robust_ci,
                                                    hd_1$orig_ci) + ggtitle("e=1") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_2 <- honest_did(es = es_u4, 
                   e = 2,
                   type="relative_magnitude")
hd_2$robust_ci <- hd_2$robust_ci[-1,]
hd2plot <- createSensitivityPlot_relativeMagnitudes(hd_2$robust_ci,
                                                    hd_2$orig_ci) + ggtitle("e=2") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_3 <- honest_did(es = es_u4, 
                   e = 3,
                   type="relative_magnitude")
hd_3$robust_ci <- hd_3$robust_ci[-1,]
hd3plot <- createSensitivityPlot_relativeMagnitudes(hd_3$robust_ci,
                                                    hd_3$orig_ci) + ggtitle("e=3") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_4 <- honest_did(es = es_u4, 
                   e = 4,
                   type="relative_magnitude")
hd_4$robust_ci <- hd_4$robust_ci[-1,]
hd4plot <- createSensitivityPlot_relativeMagnitudes(hd_4$robust_ci,
                                                    hd_4$orig_ci) + ggtitle("e=4") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_5 <- honest_did(es = es_u4, 
                   e = 5,
                   type="relative_magnitude", grid.ub = .4, grid.lb=-.4)
hd_5$robust_ci <- hd_5$robust_ci[-1,]
hd5plot <- createSensitivityPlot_relativeMagnitudes(hd_5$robust_ci,
                                                    hd_5$orig_ci) + ggtitle("e=5") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_6 <- honest_did(es = es_u4, 
                   e = 6,
                   type="relative_magnitude", grid.ub = .4, grid.lb=-.4)
hd_6$robust_ci <- hd_6$robust_ci[-1,]
hd6plot <- createSensitivityPlot_relativeMagnitudes(hd_6$robust_ci,
                                                    hd_6$orig_ci) + ggtitle("e=6") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_7 <- honest_did(es = es_u4, 
                   e = 7,
                   type="relative_magnitude",  grid.ub = .4, grid.lb=-.4)
hd_7$robust_ci <- hd_7$robust_ci[-1,]
hd7plot <- createSensitivityPlot_relativeMagnitudes(hd_7$robust_ci,
                                                    hd_7$orig_ci) + ggtitle("e=7") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_8 <- honest_did(es = es_u4, 
                   e = 8,
                   type="relative_magnitude",  grid.ub = .4, grid.lb=-.4)
hd_8$robust_ci <- hd_8$robust_ci[-1,]
hd8plot <- createSensitivityPlot_relativeMagnitudes(hd_8$robust_ci,
                                                    hd_8$orig_ci) + ggtitle("e=8") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_9 <- honest_did(es = es_u4, 
                   e = 9,
                   type="relative_magnitude", grid.ub = .4, grid.lb=-.4)
hd_9$robust_ci <- hd_9$robust_ci[-1,]
hd9plot <- createSensitivityPlot_relativeMagnitudes(hd_9$robust_ci,
                                                    hd_9$orig_ci) + ggtitle("e=9") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_10 <- honest_did(es = es_u4, 
                    e = 10,
                    type="relative_magnitude",  grid.ub = .4, grid.lb=-.4)
hd_10$robust_ci <- hd_10$robust_ci[-1,]
hd10plot <- createSensitivityPlot_relativeMagnitudes(hd_10$robust_ci,
                                                     hd_10$orig_ci) + ggtitle("e=10") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_11 <- honest_did(es = es_u4, 
                    e = 11,
                    type="relative_magnitude",grid.ub = .4, grid.lb=-.4)
hd_11$robust_ci <- hd_11$robust_ci[-1,]
hd11plot <- createSensitivityPlot_relativeMagnitudes(hd_11$robust_ci,
                                                     hd_11$orig_ci) + ggtitle("e=11") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_12 <- honest_did(es = es_u4, 
                    e = 12,
                    type="relative_magnitude", grid.ub = .4, grid.lb=-.4)
hd_12$robust_ci <- hd_12$robust_ci[-1,]
hd12plot <- createSensitivityPlot_relativeMagnitudes(hd_12$robust_ci,
                                                     hd_12$orig_ci) + ggtitle("e=12") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

combined4 <- (hd1plot + hd2plot + hd3plot) /
  (hd4plot + hd5plot + hd6plot)/  (hd7plot +
                                     hd8plot + hd9plot) / (hd10plot+ hd11plot +
                                                             hd12plot)  & theme(legend.position = "bottom")
figurea12 <- combined4 + plot_layout(guides = "collect")

ggsave("Output/figurea12.pdf", width = 8, height = 8, units="in")
print("Figure A12 complete")

## 2 pre-treatment periods

es_u2 <- aggte(att_u, type = "dynamic", na.rm = T, min_e = -4, max_e = 12)
hd_1 <- honest_did(es = es_u2, 
                   e = 1,
                   type="relative_magnitude")
hd_1$robust_ci <- hd_1$robust_ci[-1,]

hd1plot <- createSensitivityPlot_relativeMagnitudes(hd_1$robust_ci,
                                                    hd_1$orig_ci) + ggtitle("e=1") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_2 <- honest_did(es = es_u2, 
                   e = 2,
                   type="relative_magnitude")
hd_2$robust_ci <- hd_2$robust_ci[-1,]
hd2plot <- createSensitivityPlot_relativeMagnitudes(hd_2$robust_ci,
                                                    hd_2$orig_ci) + ggtitle("e=2") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_3 <- honest_did(es = es_u2, 
                   e = 3,
                   type="relative_magnitude")
hd_3$robust_ci <- hd_3$robust_ci[-1,]
hd3plot <- createSensitivityPlot_relativeMagnitudes(hd_3$robust_ci,
                                                    hd_3$orig_ci) + ggtitle("e=3") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_4 <- honest_did(es = es_u2, 
                   e = 4,
                   type="relative_magnitude")
hd_4$robust_ci <- hd_4$robust_ci[-1,]
hd4plot <- createSensitivityPlot_relativeMagnitudes(hd_4$robust_ci,
                                                    hd_4$orig_ci) + ggtitle("e=4") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_5 <- honest_did(es = es_u2, 
                   e = 5,
                   type="relative_magnitude", grid.ub = .4, grid.lb=-.4)
hd_5$robust_ci <- hd_5$robust_ci[-1,]
hd5plot <- createSensitivityPlot_relativeMagnitudes(hd_5$robust_ci,
                                                    hd_5$orig_ci) + ggtitle("e=5") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_6 <- honest_did(es = es_u2, 
                   e = 6,
                   type="relative_magnitude", grid.ub = .4, grid.lb=-.4)
hd_6$robust_ci <- hd_6$robust_ci[-1,]
hd6plot <- createSensitivityPlot_relativeMagnitudes(hd_6$robust_ci,
                                                    hd_6$orig_ci) + ggtitle("e=6") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_7 <- honest_did(es = es_u2, 
                   e = 7,
                   type="relative_magnitude",  grid.ub = .4, grid.lb=-.4)
hd_7$robust_ci <- hd_7$robust_ci[-1,]
hd7plot <- createSensitivityPlot_relativeMagnitudes(hd_7$robust_ci,
                                                    hd_7$orig_ci) + ggtitle("e=7") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_8 <- honest_did(es = es_u2, 
                   e = 8,
                   type="relative_magnitude",  grid.ub = .4, grid.lb=-.4)
hd_8$robust_ci <- hd_8$robust_ci[-1,]
hd8plot <- createSensitivityPlot_relativeMagnitudes(hd_8$robust_ci,
                                                    hd_8$orig_ci) + ggtitle("e=8") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_9 <- honest_did(es = es_u2, 
                   e = 9,
                   type="relative_magnitude", grid.ub = .4, grid.lb=-.4)
hd_9$robust_ci <- hd_9$robust_ci[-1,]
hd9plot <- createSensitivityPlot_relativeMagnitudes(hd_9$robust_ci,
                                                    hd_9$orig_ci) + ggtitle("e=9") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_10 <- honest_did(es = es_u2, 
                    e = 10,
                    type="relative_magnitude",  grid.ub = .4, grid.lb=-.4)
hd_10$robust_ci <- hd_10$robust_ci[-1,]
hd10plot <- createSensitivityPlot_relativeMagnitudes(hd_10$robust_ci,
                                                     hd_10$orig_ci) + ggtitle("e=10") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_11 <- honest_did(es = es_u2, 
                    e = 11,
                    type="relative_magnitude",grid.ub = .4, grid.lb=-.4)
hd_11$robust_ci <- hd_11$robust_ci[-1,]
hd11plot <- createSensitivityPlot_relativeMagnitudes(hd_11$robust_ci,
                                                     hd_11$orig_ci) + ggtitle("e=11") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

hd_12 <- honest_did(es = es_u2, 
                    e = 12,
                    type="relative_magnitude", grid.ub = .4, grid.lb=-.4)
hd_12$robust_ci <- hd_12$robust_ci[-1,]
hd12plot <- createSensitivityPlot_relativeMagnitudes(hd_12$robust_ci,
                                                     hd_12$orig_ci) + ggtitle("e=12") +
  theme_minimal()+ theme(legend.position="bottom") + scale_color_discrete(name="")

combined2 <- (hd1plot + hd2plot + hd3plot) /
  (hd4plot + hd5plot + hd6plot)/  (hd7plot +
                                     hd8plot + hd9plot) / (hd10plot+ hd11plot +
                                                             hd12plot)  & theme(legend.position = "bottom")
figurea13 <- combined2 + plot_layout(guides = "collect")

ggsave("Output/figurea13.pdf", width = 8, height = 8, units="in")
print("Figure A13 complete")

# Figure A14 (Alternative estimates of Reparations on Voter Registration for Surviving Victims
# Approved for Reparations. Coefficient Plot)

eventPlotresults <- staggered(df=panel,
                              i="ID",
                              t = "time.id",
                              g= "firstrep",
                              y="vote",
                              estimand="eventstudy",
                              eventTime = -12:12)
eventPlotresults$model <- "Roth and Sant'Anna (2023)"
estimate <- c(es$att.egt)
se <- c(es$se.egt)
eventTime <- c(es$egt)
eventPlotresults1 <- as.data.frame(cbind(estimate, se, eventTime))
eventPlotresults1$model <- "Callaway and Sant'Anna (2021), Varying"
eventPlotresults1$se_neyman <- NA

# Sun and Abraham
eventPlotresults2 <- staggered_sa(df = panel, 
                                  i = "ID",
                                  t = "time.id",
                                  g = "firstrep",
                                  y = "vote", 
                                  estimand = "eventstudy",
                                  eventTime=-12:12)

eventPlotresults2$model <- "Sun and Abraham (2021)"

# Universal baseline
estimate <- c(es_u$att.egt)
se <- c(es_u$se.egt)
eventTime <- c(es_u$egt)
eventPlotresults3 <- as.data.frame(cbind(estimate, se, eventTime))
eventPlotresults3$model <- "Callaway and Sant'Anna (2021), Universal"
eventPlotresults3$se_neyman <- NA
eventPlotResults <- rbind(eventPlotresults, eventPlotresults1,eventPlotresults2,eventPlotresults3)


figurea14 <- eventPlotResults %>% 
  mutate(ymin_ptwise = estimate + 1.96*se,
         ymax_ptwise = estimate - 1.96*se)%>%
  ggplot(aes(x=eventTime, y =estimate,color=model)) +
  geom_pointrange(aes(ymin = ymin_ptwise, ymax = ymax_ptwise),
                  position = position_dodge(width = 1), size=.4)+ 
  geom_hline(yintercept =0) +
  geom_vline(xintercept =0) +
  labs(x="Event Time",y="Estimate",color="") +
  ggtitle("") +
  scale_color_manual(values=c("black","purple","#69b3a2","darkgrey"))+
  theme_minimal() +
  scale_x_continuous(limits = c(-12.6,12.6),
                     breaks = seq(-12,12,1)) +
  scale_y_continuous(breaks = function(z) seq(-.1, range(z)[2], by = 0.02))

ggsave("Output/figurea14.pdf", width = 12, height = 6, units="in")
print("Figure A14 complete")

# Figure A15 ("Outcome regression and doubly robust approaches.")
att.reg <- att_gt(yname = "vote",
                  tname = "time.id",
                  idname = "ID",
                  gname = "firstrep",
                  xformla = ~age + female,
                  data = panel,
                  clustervars = "ID",
                  est_method = "reg",
                  control_group = "notyettreated") 

es.reg <- aggte(att.reg, type = "dynamic", na.rm = T, min_e = -12, max_e = 12)

estimate <- c(es.reg$att.egt)
se <- c(es.reg$se.egt)
eventTime <- c(es.reg$egt)
eventPlotresults <- as.data.frame(cbind(estimate, se, eventTime))
eventPlotresults$model <- " Outcome regression"

estimate <- c(es$att.egt)
se <- c(es$se.egt)
eventTime <- c(es$egt)
eventPlotresults1 <- as.data.frame(cbind(estimate, se, eventTime))
eventPlotresults1$model <- "Doubly robust"
eventPlotResults <- rbind(eventPlotresults, eventPlotresults1)

figurea15 <- eventPlotResults %>% 
  mutate(ymin_ptwise = estimate + 1.96*se,
         ymax_ptwise = estimate - 1.96*se)%>%
  ggplot(aes(x=eventTime, y =estimate,color=model)) +
  geom_pointrange(aes(ymin = ymin_ptwise, ymax = ymax_ptwise),
                  position = position_dodge(width = 1), size=.4)+ 
  geom_hline(yintercept =0) +
  #geom_vline(xintercept =0) +
  labs(x="Event Time",y="Estimate",color="") +
  ggtitle("") +
  scale_color_manual(values=c("black","purple"))+
  theme_minimal() +
  scale_x_continuous(limits = c(-12.6,12.6),
                     breaks = seq(-12,12,1)) +
  scale_y_continuous(breaks = function(z) seq(-.1, range(z)[2], by = 0.02))

ggsave("Output/figurea15.pdf", width = 12, height = 6, units="in")
print("Figure A15 complete")

## Figure A16 ("Estimates of Impact of Reparations on Voter Registration for Surviving
#Victims Approved for Reparations using Matching Method.")
## Matching method (panelmatch)
panel.df <- as.data.frame(panel)
panel.match <- panel.df %>%
  dplyr::select(age, female, vote, time.id, ID, firstrep, reparation)
panel.match <- panel.match[order(panel.match$ID, panel.match$time.id), ] # Order

PM.results <- PanelMatch(
  lag = 1, time.id = "time.id", unit.id = "ID",
  treatment = "reparation", refinement.method = "mahalanobis",
  qoi = "att", outcome.var = "vote",
  match.missing = F, lead=0:12,
  covs.formula = ~ age + female, forbid.treatment.reversal = TRUE,
  size.match = 1,   data = panel.match)

PE.results <- PanelEstimate(
  sets = PM.results,
  data = panel.match)
pe.results <- summary(PE.results)
coefs <- c(pe.results$summary[c(1:13)])
ses <- c(pe.results$summary[14:26])
months <- c(0:12)
coefs.df <- as.data.frame(cbind(coefs, ses, months))
figurea16 <- ggplot(coefs.df, aes(months,coefs)) + 
  geom_hline(yintercept=0, lty=2, lwd=0, colour="grey50") +
  geom_point(size=3, color = "firebrick4") +
  geom_errorbar(aes(ymin=coefs - 1.96*ses, ymax=coefs + 1.96*ses), 
                lwd=.5, width=0, color = "firebrick4") +
  theme_minimal() +
  labs(x="Months to approval",y="Change in registering") +
  scale_x_continuous(limits = c(0,12),
                     breaks = seq(0,12,1)) +
  scale_y_continuous(breaks = function(z) seq(0, range(z)[2], by = 0.01))

ggsave("Output/figurea16.pdf", width = 10, height = 5, units="in")
print("Figure A16 complete")


# Figure A17
## ATT by income bucket
lowinc.nat <- panel %>%
  filter(nat.tercile=="Low")

lowerincome.nat.att <- att_gt(yname = "vote",
                          tname = "time.id",
                          idname = "ID",
                          gname = "firstrep",
                          xformla = ~age+ female,
                          data = lowinc.nat,
                          control_group = "notyettreated")

es.lower.nat <- aggte(lowerincome.nat.att, type = "dynamic", na.rm = T, min_e = -12, max_e = 12)

midinc.nat <- panel %>%
  filter(nat.tercile=="Middle")

midincome.nat.att <- att_gt(yname = "vote",
                        tname = "time.id",
                        idname = "ID",
                        gname = "firstrep",
                        xformla = ~age+ female,
                        data = midinc.nat,
                        control_group = "notyettreated")

es.middle.nat <- aggte(midincome.nat.att, type = "dynamic", na.rm = T, min_e = -12, max_e = 12)

highinc.nat <- panel %>%
  filter(nat.tercile=="High")

highincome.att.nat <- att_gt(yname = "vote",
                         tname = "time.id",
                         idname = "ID",
                         gname = "firstrep",
                         xformla = ~age + female,
                         data = highinc.nat,
                         control_group = "notyettreated")

es.high.nat <- aggte(highincome.att.nat, type = "dynamic", na.rm=T, min_e=-12, max_e= 12)

low.nat.coefs <- make_coefs_dfs(es.lower.nat, "Low income")
mid.nat.coefs <- make_coefs_dfs(es.middle.nat, "Middle income")
hi.nat.coefs <- make_coefs_dfs(es.high.nat, "High income")

figurea17 <- make_coef_plot(low.nat.coefs, mid.nat.coefs, hi.nat.coefs) + 
  scale_y_continuous(limits = c(-.25,.22), breaks = seq(-.25, .22, by = 0.02))

ggsave("Output/figurea17.pdf", width = 10, height = 5, units="in")
print("Figure A17 complete")

# Figure A18
lowinc.nat.half <- panel %>%
  filter(nat.half=="Low")

lowerincome.nat.half.att <- att_gt(yname = "vote",
                              tname = "time.id",
                              idname = "ID",
                              gname = "firstrep",
                              xformla = ~age+ female,
                              data = lowinc.nat.half,
                              control_group = "notyettreated")

es.lower.half <- aggte(lowerincome.nat.half.att, type = "dynamic", na.rm = T, min_e = -12, max_e = 12)

highinc.nat.half <- panel %>%
  filter(nat.half=="High")

highincome.nat.half.att <- att_gt(yname = "vote",
                         tname = "time.id",
                         idname = "ID",
                         gname = "firstrep",
                         xformla = ~age + female,
                         data = highinc.nat.half,
                         control_group = "notyettreated")

es.high.half <- aggte(highincome.nat.half.att, type = "dynamic", na.rm=T, min_e=-12, max_e= 12)

low.half.coefs <- make_coefs_dfs(es.lower.half, "Low income")
hi.half.coefs <- make_coefs_dfs(es.high.half, "High income")

figurea18 <- make_coef_plot(low.half.coefs, hi.half.coefs) 

ggsave("Output/figurea18.pdf", width = 10, height = 5, units="in")
print("Figure A18 complete")

# ------------------------------------------------------------------------------
# Complete tables included in Dataverse (Robustness)
# ------------------------------------------------------------------------------

# Table S2 ("Modeling the Effects of Reparations on Voter Registration. Balanced sample, never-
# treated control group, and first year only.")

# Balanced sample
es.balanced.table <- tidy(es.balanced)
es.balanced.table <- es.balanced.table[39:63,]
balanced <- extractFromDataFrame(es.balanced.table)

# Never-treated sample
nevertreated <- extractFromDataFrame(tidy(es_n_1))
nt_att <- nevertreated[[1]]@coef
nt_se <- nevertreated[[1]]@se

# First year only
yr1 <- extractFromDataFrame(tidy(es.yr1))
yr1_att <- yr1[[1]]@coef
yr1_se <- yr1[[1]]@se

tables2 <- texreg(c(balanced, nevertreated, yr1), digits=3, custom.note= "", caption = "", label = "",
                  custom.gof.rows = list("Covariates" = c("Yes","Yes","Yes"),
                                         Individuals = c(es.balanced$DIDparams$n, es_n_1$DIDparams$n,
                                                         es.yr1$DIDparams$n),
                                         Baseline = c("Varying", "Varying", "Varying")))

event_time <- -12:12

tables2_models <- list(
  balanced = list(
    e = event_time,
    att = es.balanced.table$estimate,
    se = es.balanced.table$std.error,
    n = es.balanced$DIDparams$n
  ),
  nevertreated = list(
    e = event_time <- -12:12,
    att = as.numeric(nt_att),
    se = as.numeric(nt_se),
    n = es_n_1$DIDparams$n
  ),
  yr1 = list(
    e = event_time <- -12:12,
    att = as.numeric(yr1_att),
    se = as.numeric(yr1_se),
    n = es.yr1$DIDparams$n
  )
)

# Save as RDS
saveRDS(tables2_models, file = "Output/tables2.rds")

print("Table S2 complete")

# Table S3 ("Modeling the Effects of Reparations on Voter Registration. Longer pre-treatment period.")
prelong <- extractFromDataFrame(tidy(es.prelong))

tables3 <- texreg(prelong, digits=3, custom.note= "", caption = "", label = "",
                  custom.gof.rows = list("Covariates" = "Yes",
                                         Individuals = es.prelong$DIDparams$n,
                                         Baseline = "Universal"))

tables3_model <- list(
  prelong = list(
    e = -48:12,
    att = tidy(es.prelong)$estimate,
    se = tidy(es.prelong)$std.error,
    n = es.prelong$DIDparams$n
  )
)

saveRDS(tables3_model, file = "Output/tables3.rds")

print("Table S3 complete")

# Table S4 ("Modeling the Effects of Reparations on Voter Registration. Longer post-treatment period.")
postlong <- extractFromDataFrame(tidy(es.postlong))

tables4 <- texreg(postlong, digits=3, custom.note= "", caption = "", label = "",
                  custom.gof.rows = list("Covariates" = "Yes",
                                         Individuals = es.postlong$DIDparams$n,
                                         Baseline = "Varying"))
tables4_model <- list(
  postlong = list(
    e = -12:48,
    att = tidy(es.postlong)$estimate,
    se = tidy(es.postlong)$std.error,
    n = es.postlong$DIDparams$n
  )
)

saveRDS(tables4_model, file = "Output/tables4.rds")

print("Table S4 complete")

# Table S5 ("Modeling the Effects of Reparations on Voter Registration. Calendar time.")
calendar <- extractFromDataFrame(tidy(calendar.es))

tables5 <- texreg(calendar, digits=3, custom.note= "", caption = "", label = "",
                  custom.gof.rows = list("Covariates" = "Yes",
                                         Individuals = calendar.es$DIDparams$n,
                                         Baseline = "N/A"))

tables5_model <- list(
  calendar = list(
    e = tidy(calendar.es)$term,
    att = tidy(calendar.es)$estimate,
    se = tidy(calendar.es)$std.error,
    n = calendar.es$DIDparams$n
  )
)

saveRDS(tables5_model, file = "Output/tables5.rds")

print("Table S5 complete")

# Table S6 ("Alternative estimates of Reparations on Voter Registration for 
# Surviving Victims Approved for Reparations.")
# Roth and Sant'Anna, Sun and Abraham
rothsantanna <- staggered(df=panel,
                          i="ID",
                          t = "time.id",
                          g= "firstrep",
                          y="vote",
                          estimand="eventstudy",
                          eventTime = -12:12)

sunabraham <- staggered_sa(df = panel, 
                           i = "ID",
                           t = "time.id",
                           g = "firstrep",
                           y = "vote", 
                           estimand = "eventstudy",
                           eventTime=-12:12)

tr.rs <- createTexreg(coef = rothsantanna$estimate, coef.names = as.character(rothsantanna$eventTime), 
                      se = rothsantanna$se)
tr.sa <- createTexreg(coef = sunabraham$estimate, coef.names = as.character(sunabraham$eventTime), 
                      se = sunabraham$se)

tables6 <- texreg(c(tr.rs, tr.sa),  digits=3, custom.note= "", 
                  custom.model.names = c("Roth & Sant'Anna (2023)", "Sun and Abraham (2021)"),
                  custom.gof.rows = list("Covariates" = c("No", "No")))

tables6_model <- list(
  rothsantanna = list(
    e = rothsantanna$eventTime,
    att = rothsantanna$estimate,
    se = rothsantanna$se
  ),
  sunabraham = list(
    e = sunabraham$eventTime,
    att = sunabraham$estimate,
    se = sunabraham$se
  )
)

saveRDS(tables6_model, file = "Output/tables6.rds")

print("Table S6 complete")

## Table S7 ("Modeling the Effects of Reparations on Voter Registration.
# Outcome regression approach.")

or <- extractFromDataFrame(tidy(es.reg))
or_att <- or[[1]]@coef
or_se <- or[[1]]@se

tables7 <- texreg(or, digits=3, custom.note= "", caption = "", label = "",
                  custom.gof.rows = list("Covariates" = "Yes",
                                         Individuals = es.reg$DIDparams$n,
                                         Baseline = "N/A"))

tables7_model <- list(
  outcome_regression = list(
    e = -12:12,
    att = or_att,
    se = or_se,
    n = es.reg$DIDparams$n
  )
)

saveRDS(tables7_model, file = "Output/tables7.rds")

print("Table S7 complete")

## Table S8 ("Estimates of Impact of Reparations on Voter Registration for
# Surviving Victims Approved for Reparations using Matching Method.")
coefs <- c(pe.results$summary[c(1:13)])
ses <- c(pe.results$summary[14:26])
months <- c(0:12)
coefs.df <- as.data.frame(cbind(coefs, ses, months))

pm <- createTexreg(coef = coefs.df$coefs, coef.names = as.character(coefs.df$months), 
                   se = coefs.df$ses)

tables8 <- texreg(pm, digits=3, custom.note= "", caption = "", label = "",
                  custom.gof.rows = list("Covariates" = "Yes"))

tables8_model <- list(
  matching_method = list(
    e = coefs.df$months,
    att = coefs.df$coefs,
    se = coefs.df$ses
  )
)

saveRDS(tables8_model, file = "Output/tables8.rds")

print("Table S8 complete")

## Table S9 ("Effects of Reparations on Voter Registration by Municipality
# Income Level (defined by national tercile).")
col1 <- extractFromDataFrame(tidy(es.lower.nat))
col2 <- extractFromDataFrame(tidy(es.middle.nat))
col3 <- extractFromDataFrame(tidy(es.high.nat))

tables9 <- texreg(c(col1, col2, col3),  digits=3,
                  custom.model.names = c("Low", "Middle", "High"),
                  custom.note= "", caption = "", label = "",
                  custom.header = list("Income level" = 1:3),
                  custom.gof.rows = list("Covariates" = c("Yes", "Yes", "Yes"),
                                         Individuals = c(es.lower.nat$DIDparams$n, es.middle.nat$DIDparams$n,
                                                         es.high.nat$DIDparams$n),
                                         Baseline = c("Varying", "Varying", "Varying")))

tables9_models <- list(
  low_income_national = list(
    e = es.lower.nat$egt,
    att = es.lower.nat$att.egt,
    se = es.lower.nat$se.egt,
    n = es.lower.nat$DIDparams$n
  ),
  middle_income_national = list(
    e = es.middle.nat$egt,
    att = es.middle.nat$att.egt,
    se = es.middle.nat$se.egt,
    n = es.middle.nat$DIDparams$n
  ),
  high_income_national = list(
    e = es.high.nat$egt,
    att = es.high.nat$att.egt,
    se = es.high.nat$se.egt,
    n = es.high.nat$DIDparams$n
  )
)

saveRDS(tables9_models, file = "Output/tables9.rds")

print("Table S9 complete")

## Table S10 ("Effects of Reparations on Voter Registration by Municipality
# Income Level (defined by national median).")
col1 <- extractFromDataFrame(tidy(es.lower.half))
col2 <- extractFromDataFrame(tidy(es.high.half))

tables10 <- texreg(c(col1, col2),  digits=3,
                  custom.model.names = c("Low", "High"),
                  custom.note= "", caption = "", label = "",
                  custom.header = list("Income level" = 1:2),
                  custom.gof.rows = list("Covariates" = c("Yes", "Yes"),
                                         Individuals = c(es.lower.half$DIDparams$n,
                                                         es.high.half$DIDparams$n),
                                         Baseline = c("Varying", "Varying")))

tables10_models <- list(
  low_income_half = list(
    e = es.lower.half$egt,
    att = es.lower.half$att.egt,
    se = es.lower.half$se.egt,
    n = es.lower.half$DIDparams$n
  ),
  high_income_half = list(
    e = es.high.half$egt,
    att = es.high.half$att.egt,
    se = es.high.half$se.egt,
    n = es.high.half$DIDparams$n
  )
)

saveRDS(tables10_models, file = "Output/tables10.rds")

print("Table S10 complete")

# Table S11: ("Effects of Reparations on Voter Registration by Age at Approval.")
panel.aap <- panel %>%
  mutate(age.at.first = (as.yearmon(StartDate) - as.yearmon(F_Nac)))

over70 <- panel.aap %>%
  filter(age.at.first > 70) # greater than 70

over70.model <- att_gt(yname = "vote",
                       tname = "time.id",
                       idname = "ID",
                       gname = "firstrep",
                       xformla = ~age + female,
                       data = over70,
                       control_group = "notyettreated")

es.over70 <- aggte(over70.model, type = "dynamic", na.rm = T, min_e = -12, max_e = 12)


under70 <- panel.aap %>%
  filter(age.at.first <= 70) # younger than 70

under70.model <- att_gt(yname = "vote",
                        tname = "time.id",
                        idname = "ID",
                        gname = "firstrep",
                        xformla = ~age + female,
                        data = under70,
                        control_group = "notyettreated")

es.under70<- aggte(under70.model, type = "dynamic", na.rm = T, min_e = -12, max_e = 12)

col1 <- extractFromDataFrame(tidy(es.over70))
col2 <- extractFromDataFrame(tidy(es.under70))

tables11 <- texreg(c(col1, col2),  digits=3,
                   custom.model.names = c("Over 70", "70 or below"),
                   custom.note= "", caption = "", label = "",
                   custom.gof.rows = list("Covariates" = c("Yes", "Yes"),
                                          Individuals = c(es.over70$DIDparams$n, es.under70$DIDparams$n),
                                          Baseline = c("Varying", "Varying")))
tables11_models <- list(
  over70 = list(
    e = es.over70$egt,
    att = es.over70$att.egt,
    se = es.over70$se.egt,
    n = es.over70$DIDparams$n
  ),
  under70 = list(
    e = es.under70$egt,
    att = es.under70$att.egt,
    se = es.under70$se.egt,
    n = es.under70$DIDparams$n
  )
)

saveRDS(tables11_models, file = "Output/tables11.rds")

print("Table S11 complete")

## Stat used in body of paper
# Goodman-Bacon decomposition - stat is used in the body of the paper.
# Estimated runtime: 5 hours
bacondecomp <- bacon(vote ~ reparation,
                     data = panel,
                     id_var = "ID",
                     time_var = "time.id")